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SOLVING CONSTRAINED QUADRATIC BINARY PROBLEMS 
VIA QUANTUM ADIABATIC EVOLUTION 

POOYA RONAGH, BRAD WOODS, AND EHSAN IRANMANESH 


Abstract. Quantum adiabatic evolution is perceived as useful for binary quadratic pro¬ 
gramming problems that are a priori unconstrained. For constrained problems, it is a com¬ 
mon practice to relax linear equality constraints as penalty terms in the objective function. 
However, there has not yet been proposed a method for efficiently dealing with inequality 
constraints using the quantum adiabatic approach. In this paper, we give a method for solv¬ 
ing the Lagrangian dual of a binary quadratic programming (BQP) problem in the presence 
of inequality constraints and employ this procedure within a branch-and-bound framework 
for constrained BQP (CBQP) problems. 


1. Introduction 

The unconstrained binary quadratic programming (UBQP) problem is defined by 

Minimize x^Qx 
subject to a; G {0,1}"', 

where, without loss of generality, Q G Recent advancements in quantum computing 

technology [61E01I22] have raised hopes of the production of computing systems that are 
capable of solving UBQP problems, showing quantum speedup. The stochastic nature of 
such systems, together with all sources of noise and error, are challenges yet to be overcome 
in achieving scalable quantum computing systems of this type. This paper is regardless 
motivated by the assumption of the existence of systems that can solve UBQP problems 
efficiently and to optimality, or at least in conjunction with a framework of noise analysis of 
the suboptimal results. We call such a computing system a UBQP oracle. 
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Many NP-hard combinatorial optimization problems arise naturally or can easily be re¬ 
formulated as UBQP problems, such as the quadratic assignment problem, the maximum 
cut problem, the maximum clique problem, the set packing problem, and the graph colour¬ 
ing problem (see, for instance. Boros and Prekopa [H], Boros and Hammer [S], Bourjolly et 
al. HU], Du and Pardalos [H], Pardalos and Rodgers [2ZII2B], Pardalos and Xue [2U], and 
Kochenberger et al. 1231 ). 

Numerous interesting applications that are expressed naturally in the form of UBQP 
problems appear in the literature. Barahona et al. mm formulate and solve the problem 
of hnding exact ground states of spin glasses with magnetic helds. Alidaee et al. |2] study 
the problem of scheduling n jobs non-preemptively on two parallel identical processors to 
minimize weighted mean flow time as a UBQP problem. Bomze et al. give a comprehensive 
discussion of the maximum clique (MC) problem. Included is the UBQP representation of 
the MC problem and a variety of applications from different domains. UBQP has been used 
in the prediction of epileptic seizures [TU]. Alidaee et al. PQ discuss a number partitioning 
problem, formulating a special case as a UBQP problem. 

In the literature, problems which naturally occur as constrained binary quadratic pro¬ 
gramming (CBQP) problems are commonly reformulated as UBQP problems by including 
quadratic penalties in the objective function as an alternative to explicitly imposing con¬ 
straints. Although this method has been used very successfully on classical hardware, it is 
not a viable approach when using quantum adiabatic hardware, as the reformulation dra¬ 
matically increases the density, range of coefficients, and dimension of the problem. 

In this paper, we consider the CBQP problem with linear constraints, stated formally as 

(P) Minimize x^Qx 
subject to Ax < b, 

xe{0,ir, 

where Q G and A G We present a branch-and-bound approach which uses 

Lagrangian duality to solve (P) and show that a UBQP oracle can be used to solve the 
Lagrangian dual (LD) problem with successive applications of linear programming (LP). We 
introduce the notion of quantum annealing leniency, which represents the maximum thresh¬ 
old of the annealing time allowed for our approach to outperform the Gurobi Optimizer. 

The remainder of this paper is organized as follows. Section [2] gives an overview of the 
quantum adiabatic approach to computation and the relationship of the approach to UBQP. 
Section [3] presents lower-bounding procedures for CBQP. Section 0] presents a local search 
heuristic for CBQP. Branching strategies are described in Section [5l Test instances and 
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results of our computational experiments are presented in Section [6l In Section [7] we provide 
practical remarks for application of our method to the D-Wave systems, and in Section [8] we 
mention how our method can more generally be used to solve constrained binary program¬ 
ming problems with higher-order objective and constraint polynomials. 


2. Preliminaries 

2.1. Overview of quantum computing using adiabatic evolution. Recent advance¬ 
ments in quantum hardware have increased the motivation to study forms of computation 
that differ in computational complexity from the limitations of Turing machines. The quan¬ 
tum gate model is a means of achieving powerful quantum algorithms such as Shor’s well 
known quantum algorithms for factoring and discrete logarithms |32]. Aside from the quan¬ 
tum gate model, there are several other paradigms of quantum information technology, each 
of which would open a new world of possible algorithm designs for the realization of a cor¬ 
responding practical quantum processor. 

Farhi et ah [ig M\ propose quantum adiabatic evolution as a novel paradigm for the 
design of quantum algorithms. The physicist’s interpretation of this quantum phenomenon 
is a “quantum” analogue to what an optimization theorist would view as simulated annealing. 
In particular, quantum adiabatic evolution can be viewed as a quantum local search that 
converges to a best known solution as a parameter t varies from an initial time 0 to a hnal 
time T. Just as in the case of a simulated annealing algorithm, the challenge is to show that 
the process converges to the best solution with no n zero probability if T is a polynomial in 
the size of the problem. 

Van Dam et ah [33] give an example of an adiabatic quantum algorithm for searching that 
matches the optimal quadratic speedup obtained by Grover’s search algorithm. This example 
demonstrates that the “quantum local search”, which is implicit in the adiabatic evolution, 
is truly non-classical in nature from a computational perspective. In the same reference |33j . 
Theorem 4.1 explains how the continuous-time evolution of t G [0,T] can be approximated 
by a quantum circuit consisting of a sequence of poly(?7,T) unitary transformations. 

All of the above considerations motivate the main assumption of this paper: practical 
quantum hardware can result in a signihcant quantum speedup in certain forms of integer 
programming. Our goal is to design and suggest optimization algorithms that work in 
conjunction with such integer programming oracles. 

In the rest of this section, following [33], we review more details of the quantum adiabatic 
evolution and explain the quantum adiabatic approach to UBQP. 
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2.2. The quantum adiabatic theorem. A time-dependent quantum mechanical system 
is described by the Schrodinger equation 

= H{s)\i;{s)). 

Here h is Planck’s constant, s is a reparametrization of time ranging in [0,1]- is a 

vector in a high-dimensional Hilbert space and represents the state of the system. H{s) 
is the Hamiltonian of the system, a Hermitian operator on the Hilbert space such that its 
eigenvectors form the set of all possible states of the system, also known as eigenstates. The 
eigenvalue spectrum of H{s) refer to the different energies of the eigenstates. The state 
(eigenvector) with the lowest energy (eigenvalue) is called the ground state of the system. 

In most scenarios the Schrodinger equation is a complicated differential equation to tackle. 
The main challenge in quantum physics is approximating and understanding the properties 
of the solution of this equation. The quantum adiabatic theorem is an example of such 
investigation. 

The theorem states that a physical system that is initially in its ground state tends to stay 
in that state, provided that the Hamiltonian of the system is changed slowly enough |2B] . 
For instance, if we add a constant delay factor T to the Schrodinger equation. 

Minis)) = TH(s)Ws)), 

the evolution is slow enough in the sense of the adiabatic theorem, provided that T satishes 

T > , for all se0,l. 

9{sy 

Here g{s) is the difference between the smallest two eigenvalues of H{s) and is often referred 
to as the gap at time s. Finding g as a function of s is a difficult problem of interest in 
condensed matter physics. However, to give a constant delay factor we need only to hnd a 
lower bound for it. We may therefore simply consider the minimum gap := min,, g{s) and 
also Amax := rnaxs ||^i/(s)|| 2 . We can then verify that a constant delay factor T G 0{^^) 

drain 

is sufficient for the adiabatic evolution. 


2.3. Adiabatic quantum computation. Let / : {0,1}"^ —)■ M be a function on the binary 
domain that is computable in polynomial time and bounded in values by a polynomial in 
the size of its input argument (for example, / is a polynomial). In adiabatic quantum 
computation as proposed by Farhi et ah [15], we associate a physical system to / such that 
its Hamiltonian Hf takes the same eigenvalue spectrum as /. Hence Hj can be explicitly 
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written as 


Hf ■■= Y1 

ze{o,i}" 

in the basis of its eigenstates, also encoded as the domain {0,1}”. 

We consider another Hamiltonian Ho whose eigenstates are easy to compute but is not 
diagonal in the z-basis. An explicit suggestion for Hq is as follows. Given 

|0):=;^(|0) + |1)) and |i) ;= ^(|0) - |1)), 

we consider the “Hadamard basis” as consisting of all states \z) which would be written as 
\z) in this basis as z ranges over all binary strings z G {0,1}". Then Hq may be dehned as 

Ho := ^ h{z)-\z){z\, 

2e{o,i}" 

with = 0 and h{z) > 1 for all other z ^ 0”. Thus the ground state of Ho is |0 ■ ■ -0), 
which is a uniform superposition of all of the states in the z-basis, 

Z 

For any practical implementation it is necessary that the hardware can be easily set to this 
ground state. 

The quantum adiabatic computation is dehned as the time-dependent Hamiltonian 

Hit) ;= {l-L)Ho + ^Hf, 

with 0 < f < T when the system is set to the ground state of Ho at time t = 0: 

1^(0)) = |o-). 

Here T > 0 is the delay factor of the evolution and by the adiabatic theorem this system 
will evolve |'0(O)) to the global minimum of the function /, provided that T G 

^min 

In what follows we work under the assumption of the existence of a UBQP oracle, an 
oracle for solving UBQP problems. This is specihcally motivated by current progress of 
quantum adiabatic technology pioneered by D-Wave Systems, in which the connectivities of 
the quantum bits are only in the form of couplings between pairs of quantum bits [20] . How¬ 
ever, we will observe that our suggested methods are easily generalizable to take advantage 
of systems with higher-degree interactions of quantum bits if such systems are implemented 
in the future (see Section [H|). 
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3. Lower-bounding procedures 


3.1. Linearization relaxation. A standard linearization of (P) involves relaxing the in¬ 
tegrality constraint on variables Xj, (z = 1,... ,n) and defining continnous variables yij for 
every pair XiXj in the objective fnnction with i < j, yielding the following linearized problem. 

n 

(Plp) Minimize 

l<2<j'<n 2=1 

snbject to Ax < b, 

Vij > Xi + Xj — 1 (V i,j snch that i < j and qij > 0), 

Vij < Xi (yi,j snch that i < j and Qij < 0), 

Vij < Xj (Vz, i snch that i < j and q^j < 0), 

0 < Xj < 1 (z = 1,..., n), 

2 /> 0 . 

A lower bonnd to (P) can now be obtained by solving (Plp) using linear programming. We 
employed this linearization in our computational experiments (see Section |6]). We recall 
that there are several methods for linearization of (P) to LP problems, many of which have 
been mentioned in the survey im by Floudas and Gounaris, and depending on the specihc 
application at hand, may be better choices. We will ignore this variety of choices, and call 
(Plp) the LP relaxation of (P). 


3.2. Lagrangian dual. We can give a lower bound for (P) by the LD problem 

(L) max d(A), 

AeR!p 

where d{\) is evaluated via the Lagrangian relaxation 

(La): d(X)= min L(x, A) = x'^Qx + X^{Ax — b). 

a:e{0,l}" 

The function d{X) is the minimum of a hnite set of linear functions of A and hence it is 
concave and piecewise linear. 

A number of techniques to solve (L) exist in the literature; however, hnding this bound is 
computationally expensive, so looser bounds (for example, the LP relaxation) are typically 
used. We note that the problem yields a natural solution using the UBQP oracle via the 
outer Lagrangian linearization method. The book by Li and Sun [25] provides background 
and several details of this approach. 
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We recall that (L) can be rewritten as an LP problem: 

(Llp) Maximize ^ 

subject to // < x^Qx + {Ax — b) (Vx G {0,1}”), 

A > 0. 

This formulation is difficult to solve directly, as there are an exponential number of con¬ 
straints. In particular, there is one linear constraint (cutting plane) for every binary point 
X G {0,1}”. However, the restriction of the constraints to a much smaller nonempty subset 
T C {0,1}” of binary vectors is a tractable LP problem: 

(Llp-t) Maximize /i 

subject to /i < x^Qx + {Ax — h) (Vx G T), 

A > 0. 

This LP problem is bounded provided that T contains at least one feasible solution. In several 
applications hnding at least one feasible solution is easy but sometimes very difficult. In the 
latter cases we impose an upper bound on the vector of Lagrange multipliers, substituting the 
constraint A>0by0<A<ri. The vector u of upper bounds, may depend on an estimation 
of the solution to (L). In practical situations, it should also depend on the specihcs of the 
UBQP oracle (for example, the noise and precision of the oracle). Let (/i*. A*) be an optimal 
solution to (Llp_t)- It is clear that (L^*) is a UBQP problem that can be solved using the 
UBQP oracle. By successively adding the solutions returned by the UBQP oracle as cutting 
planes, and a successive application of the simplex method, we are able to solve (Llp) in a 
similar fashion to [251 Procedure 3.2]. 

4. A LOCAL SEARCH HEURISTIC 

In order to prune a branch-and-bound tree effectively, it is important to quickly obtain 
a good upper bound. We employ an adaptation of the local search algorithm presented by 
Bertsimas et ah [5]. The main idea is as follows: beginning with a feasible solution x, we 
iteratively improve the solution by considering solutions in the 1-flip neighbourhood of x 
(dehned as the set of solutions that can be obtained by flipping a single element of x) which 
are feasible, together with “interesting” solutions, until it cannot be improved further. The 
algorithm takes a parameter p which explicitly controls the trade-off between complexity 
and performance by increasing the size of the neighbourhood. A neighbouring solution y is 
considered “interesting” if it satishes the following conditions: (i) no constraint is violated 
by more than one unit, and (ii) the number of violated constraints in y plus the number of 
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loose constraints which differ from the loose constraints in the cnrrent best solntion is at 
most p. 


Algorithm 1 Local search henristic 

Require: matrix A, Q; vector 6; feasible solntion zo] scalar parameter p > 0 
Ensure: Feasible solution z such that z'^Qz < ZqQzq 
1: 2 ; := Zq; S := {z} 

2 : while S' 7 ^ 0 do 
3: get a new solution x from S 

4: for all y adjacent to x do 

5: if y is feasible and y^Qy < z^Qz then 

6: z y 

7: S^y 

8: goto Step 3 

9: else if y is interesting then 

10: S^SU{y} 

11 : end if 

12 : end for 

13: end while 


We note that the algorithm moves to a better solution as soon as it finds a feasible one, 
and only when none are found does it consider moving to “interesting” solutions. 


5. Branching strategies 

In any branch-and-bound scheme, the performance of the algorithm is largely dependent 
on the number of nodes that are visited in the search tree. As such, it is important to 
make effective branching decisions, reducing the size of the search tree. Branching heuristics 
are usually classihed as either static variable-ordering (SVO) heuristics or dynamic variable¬ 
ordering (DVO) heuristics. All branching heuristics used in this paper are DVO heuristics, as 
they are generally considered more effective because they allow information obtained during 
a search to be utilized to guide the search. 

It is often quite difficult to find an assignment of values to variables that satisfies all 
constraints. This has motivated the study of a variety of approaches that attempt to ex¬ 
ploit the interplay between variable-value assignments and constraints. Examples include 
the impact-based heuristics proposed by Refalo the conflict-driven variable ordering 
heuristic proposed by Boussemart et ah and the approximated counting-based heuris¬ 
tics proposed by Kask et ah [21], Hsu et ah [18], Bras et ah [21], and Pesant et ah [30] . 






5.1. Counting the solution density. One branching heuristic used in this paper is a 
modihed implementation of the maxSD heuristic introduced by Pesant et ah [30]. We recall 
two dehnitions. 

Definition 5.1. Given a constraint c{xi,... ,x„) and respective finite domains Di [1 < i < 
n), let fi=c{xi ,..., Xn) denote the number of n-tuples in the corresponding relation. 

Definition 5.2. Given a constraint c(a;i,... ,Xn), respective finite domains Di [1 < i < n), 
a variable x* in the scope of c, and a value d G D^, the solution density of a pair (xj, d) in c 
is given by 

^c(xi,..., Xj_i, d, Xj-|_i,..., x,^) 

#c(xi, . . . ,Xn 

The solution density measures how frequently a certain assignment of a value in a variable’s 
domain belongs to a solution that satishes constraint c. 

The heuristic maxSD simply iterates over all of the variable-value pairs and chooses the 
pair that has the highest solution density. If the (approximate) a{xi,d,c) are precomputed, 
the complexity of the modihed algorithm is 0{mq), where m is the number of constraints, 
and q is the sum of the number of variables that appear in each constraint. Pesant et al. [30] 
detail good approximations of the solution densities for knapsack constraints, which can be 
computed efficiently. They also provide an in-depth experimental analysis that shows that 
this heuristic is state of the art among counting-based heuristics. 



a{xi,d, c) 


5.2. Constraint satisfaction via an LD solution. In each node u of the branch-and- 
bound tree, a lower bound is computed by solving the LD problem, (L), and the primal-dual 
pair (x“, A“) is obtained. In the standard way, we dehne the slack of constraint z at a point 
X to be Si = bi — ajx, where a* is the zth row of A. Then the set of violated constraints at x 
is the set V = {i \ Si < 0}. If x“ is infeasible for the original problem, it must violate one or 
more constraints. Additionally, we dehne the change in slack for constraint i resulting from 
hipping variable j in x“ to be 

5ij = aij{2x'f - 1) . 

We present three branching strategies which use this information at x“ to guide variable and 
value selection towards feasibility. 

The hrst branching method we propose is to select the variable that maximizes the reduc¬ 
tion in violation of the most violated constraint. That is, we select j = argmaxhjj and value 

1 — Xj (see Algorithm [2]) . 
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Algorithm 2 Most violated constraint satisfaction branching scheme 

Require: is the optimal solntion to (L) at the cnrrent node. 

1: for all constraints i do 
2: compnte Si = bi — ajx^ 

3: end for 
4: i = arg min s* 

5: if Sj > 0 then 

6: LD optimal is feasible; abort violation branching 

7: end if 

8: for all variables j do 
9: compnte 6ij = aij(2x“ — 1) 

10 : end for 

11: return the index j = argmaxhjj and valne 1 — 


The next branching method we discnss is more general: instead of looking only at the 
most violated constraint, we consider all of the violated constraints and select the variable 
which, when flipped in the LD solntion, gives the maximnm decrease in the left-hand side 
of all violated constraints (see Algorithm [3]). 


Algorithm 3 All violated constraints satisfaction branching scheme 

Require: is the optimal solution to (L) at the current node. 

1: for all constraints i do 
2: compute Si = bi — ajx^ 

3: end for 

4: 1/ := {z : Sj < 0} defines the set of violated constraints. 

5: if D = 0 then 

6: LD optimal is feasible; abort violation branching 

7: end if 

8: for all variables j do 
9: compute Sij = aij{2x'^ — 1) 

10: end for 

11: return the index j = arg max %) value 1 — 


This can be generalized to the entire constraint matrix such that we are looking for a 
variable which, when flipped in the LD solution, has the minimum increase in the left-hand 
side of the non-violated constraints and the maximum decrease in the left-hand side of the 
violated constraints (see Algorithm 0]) . 
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Algorithm 4 All constraints satisfaction branching scheme 

Require: is the optimal solntion to (L) at the cnrrent node. 

1: for all constraints i do 
2: compnte Si = bi — ajx^ 

3: end for 

4: for all constraints i do 
5: for all variables j do 

6: compnte Sij = — 1) 

7: end for 

8 : end for 



9: return 


5.3. LP-based branching. Introduced in CPLEX 7.5 [3], the idea of strong branching is 


to test which fractional variable gives the best bound before branching. The test is done 
by temporarily fixing each fractional variable to 0 or 1 and solving the LP relaxation by 
the dual simplex method. Since the cost of solving several LP subproblems is high, only a 
hxed number of iterations of the dual simplex algorithm are performed. The variable that 
provides the strongest bound is chosen as the branching variable. 

If the number of fractional variables is large, this process is very time consuming. There 
are several possible methods to overcome this difficulty. First, we can perform this test on 
only a subset of fractional variables (finding a proper subset is another issue; one possible 
way is to look only for variables with values close to 0.5). Another approach, the k-look-ahead 
branching strategy, requires an integer parameter k and a score assignment on the fractional 
variables (or variable-value pairs). The variables (or variable-value pairs) are then sorted 
according to their scores, and the above test is performed. If no better bound is found for k 
successive variables, the test process is stopped. 

We note that the lower-bound computation performed by the UBQP oracle is fully paral- 
lelizable with strong branching or the /c-look-ahead strategy performed on a digital processor, 
affording the computational time required to utilize this type of branching without signih- 
cantly increasing the total running time of the algorithm. 

5.4. Frequency-based branching. Motivated by the notion of persistencies as described 
by Boros and Hammer |H], and observing that the outer Lagrangian linearization method 
yields a number of high-quality solutions, one can perform /c-look-ahead branching, selecting 
variable-value pairs based on their frequency. Here, given a set S' C {0,1}" of binary vectors, 
an index i G {1,..., n}, and a binary value s G {0,1}, the frequency of the pair (x,, s) in S 
is defined as the number of elements in S with their ith entry equal to s. Moreover, when 
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using a UBQP oracle that performs quantum annealing, the oracle returns a spectrum of 
solutions and all solutions can be used in the frequency calculation. This branching strategy 
has not, to our knowledge, previously appeared in the literature. 

6. Computational experiments 

6.1. Generation of test instances. In this paper we use randomly generated test instances 
with inequality constraints. For a desired test instance, we let n denote the number of 
variables and m be the number of inequality constraints. To generate the cost matrix, we 
construct an n x n random symmetric matrix Q of a given density d. An m x n constraint 
matrix is generated in a similar manner, ensuring that the CBQP problem has at least one 
feasible solution. For our test instances, the numbers of constraints were chosen as half the 
numbers of variables, and densities 0.3 and 0.5 were used for the objective functions and the 
constraints matrices, respectively. 

6.2. Computational results. We now present the details of our computational experi¬ 
ments. The algorithms were programmed in C-|--|- and compiled using GNU GCC on a 
machine with a 2.5 GHz Intel Gore i5-3210M processor and 16 GB of RAM. Linear program¬ 
ming and UBQP problems were solved by the Gurobi Optimizer 6.0. We used the Gurobi 
Optimizer in place of a UBQP oracle and replaced the computational time with 0 millisec¬ 
onds per solver call. The algorithm was coded to utilize 4 cores and to allow us to accurately 
report times. All other processes were paused during the solution of the UBQP problems. 

In Tables [T] and [2] we report computational experiments performed on the group of test 
instances, evaluating each of the different branching strategies. In these tables the columns 
mostviol, allviol, and allcst respectively correspond to Algorithms [2], [3l and 01 lp4 and 
lp8 correspond respectively to the LP-based 4-look-ahead and 8-look-ahead strategies, and 
freq4 and freqS correspond respectively to the frequency-based 4-look-ahead and 8-look¬ 
ahead strategies. Table [1] gives the number of nodes that the branch and bound requires 
when using each of the branching strategies. The hnal column reports the number of nodes 
that the Gurobi Optimizer used when solving the problem directly. In terms of the number 
of nodes explored in the branch-and-bound tree, the most violated (Algorithm |2]) and all 
violated (Algorithm [3]) constraints satisfaction branching schemes are clear winners. 

Table [2] reports the time taken, in seconds, for each of the branching strategies and the 
Gurobi Optimizer to solve the problem to optimality, each using 4 cores. The number of 
queries to the UBQP oracle is given, as well as the quantum annealing leniency (QAL). 
The QAL column is computed by taking the difference between the time taken by the best 
branching strategy and the Gurobi Optimizer, and dividing by the number of queries. This 
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column 

perform 


can be viewed as the maximum threshold of the average time, in milliseconds, to 
each quantum annealing process in order to solve the original problem faster than 

the Gurobi Optimizer. We note that the 
frequency-based branching heuristics termi¬ 
nate in the least amount of computational 
. time. 

In Figure [H we graph the values from 
the QAL column of Table [2] versus problem 
size. The dotted line plots the mean val¬ 
ues. In Figure |2l we graph the values from 
the QAL column of Table [2] versus problem 
. size using a logarithmic scale. These graphs 

. suggest an exponential growth in QAL with 

respect to the problem size. We interpret 
this as an indication that the difference be- 
• tween the computational time required for 

our CBQP approach (in conjuction with a 
scalable UBQP oracle) and the Gurobi Op- 
• ‘ /' timizer grows exponentially with the size of 

/' the problem. 

• / 

• ^ • 

. , 7. Discussion 

’ ' I ^ In this section we consider the specifics of 

f'\ • ’ the D-Wave systems as a physical invention 

of a UBQP oracle. We imagine any imple¬ 
mentation of quantum adiabatic computing 
would have similar limitations, so we expect 
^^^— our algorithms to be beneficial in overcom¬ 

ing them. 

Quantum adiabatic devices have not thus 
far allowed for fully connected systems of 
quantum bits. Due to this sparsity in the 
manufactured chips, the use of such quan¬ 
tum computers requires solving a minor-embedding problem, described in the following sec¬ 
tion, prior to programming the chip according to appropriate couplings and local helds 


Number of variables 


Figure 1. Quantum Anneal¬ 
ing Leniency versus Problem 
Size 
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7.1. Efficient embedding. Given a UBQP 
instance defined by a matrix Q, we define 
the underlying graph H as follows: for each 
variable Xi, we associate a vertex Vi eV (i7); 
and for each nonzero entry qij ^ 0 of Q with 
i ^ j, we let Vi and Vj be adjacent in H. 

The minor-embedding problem that must be 
solved is to find a function cj) : V{H) —)■ 

where G is the graph defined by the 
quantum chip (that is, vertices correspond 
to the quantum bits and edges correspond 
to the couplings between them) such that 

(i) for each x G V{H), the subgraph in¬ 
duced by (j){x) in G is connected; 

(ii) (f){x) and (f){y) are disjoint for all x ^ y 
in V{H)-, and 

(hi) if x and y are adjacent in H, then there 
is at least one edge between 0(x) and 

<P{y) in G. 

We note that for any induced subgraph of 
H, a minor embedding can be found simply 
by restricting (j) to the subgraph. In our CBQP approach, the constraints contribute only 
to linear terms in the Lagrangian relaxation of the problem, and hence the embedding of 
a UBQP problem at any node in the branch-and-bound tree can be found from the parent 
node by restriction. That is, our method requires solving the minor-embedding problem only 
once, at the root node of the branch-and-bound tree. 

7.2. Efficient programming of the quantum chips. In every node of the branch-and- 
bound tree, all Lagrangian relaxations generated have identical quadratic terms and only 
differ from each other in linear terms. This suggests that if reprogramming the quantum 
chip can allow for fast updates of previous setups, then the runtime of UBQP oracle queries 
can also be minimized. 

7.3. Post-processing and error analysis. Currently used quantum bits have significant 
noise. For arbitrary choices of initial and final Hamiltonians, the eigenvalues in the energy 
spectrum of the evolving Hamiltonian of the system may experience gap closures. Further¬ 
more, the measurement process of the solutions of the quantum adiabatic evolution has a 



38 40 42 44 46 48 

Number of variables 


Figure 2. Quantum Anneal¬ 
ing Leniency versus Problem 
Size (log scale) 
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stochastic nature. Each of these indicate that the solutions read from the quantum system 
are often very noisy and even after several repetitions of the process there is no guarantee 
of optimality for the corresponding UBQP problem. In order to make our method practical, 
with a proof (or at least certainty) of optimality, it is necessary to develop a framework for 
error analysis for the quantum annealer. We note that for our purposes the solution errors 
can only propagate to hnal answers in the branch-and-bound tree if the proposed lower 
bound i obtained at a node is greater than the actual lower bound £ — e and the best known 
upper bound u satishes u < i. If this situation occurs, then the proposed method incorrectly 
prunes the subtree rooted at this node. 

The quantum computing community often considers post-processing to be a procedure 
of descending from the best solution provided by the UBQP oracle to a better solution. 
However, for our purposes, the post-processing needed is only in the value of the optimal 
solution (that is, given some level of required certainty, finding an estimate of e in the 
notation above). 


8. Extension to quadratically constrained problems 

It is straightforward to generalize the method proposed here to quadratically constrained 
quadratic programming (QCQP) problems in binary variables. In fact, the Lagrangian 
relaxations of QCQP problems are also UBQP problems. The minor-embedding problem to 
be solved at the root node takes the underlying graph H as follows: for each variable a:*, we 
associate a vertex Vi E V (if); and for any pair of distinct indices i j, we let (n*, Vj) G E{H) 
if and only if the term XiXj appears with a nonzero coefficient in the objective function or 
in any of the quadratic constraints. 

Assuming future quantum annealing hardware will allow higher-degree interactions be¬ 
tween quantum bits, more general polynomially constrained binary programming problems 
could also be solved using a similar approach via Lagrangian duality. 

9. Conclusions 

Motivated by recent advancements in quantum computing technology, we have provided 
a method to tackle constrained binary programming problems using these technologies. Our 
method is a branch-and-bound algorithm where the lower bounds are computed using La¬ 
grangian duality and queries to an oracle that solves unconstrained binary programming 
problems. 

The conventional branching heuristics for integer programming problems rely on fractional 
solutions of continuous relaxations of the subproblems at the nodes of the branch-and-bound 
tree. Our lower-bounding methods are not based on continuous relaxations of the binary 
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variables. In particular, the optimal solutions of the dual problems solved at the nodes of 
the branch-and-bound tree are not fractional. We have therefore proposed several branching 
strategies that rely on integer solutions and compared their performance both in time and 
number of nodes visited, to each other and to the conventional branching strategies that rely 
on fractional solutions. 

To understand how powerful quantum computing hardware needs to be (in terms of the 
time taken for each query), we introduced the notion of quantum annealing leniency. This is, 
roughly speaking, the average time a query to the UBQP oracle would be able to spend so as 
to remain as fast as a classical algorithm for solving CBQP problems. In our computational 
experiments this classical algorithm is that of the Gurobi Optimizer. 

Finally, although we focused on quadratic objective functions and linear inequality con¬ 
straints, in Section E] we discuss how this method can be generalized to higher-order binary 
polynomial programming problems. 
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Table 1. Branch-and-bound nodes needed to solve each instance 


size 

mostviol 

allviol 

allcst 

lp4 

lp8 

freq4 

freqS 

gurobi 

38 

9 

9 

95 

115 

79 

55 

33 

14232 


19 

19 

193 

87 

35 

63 

55 

2823 


29 

33 

57 

93 

77 

267 

141 

7020 


9 

9 

69 

23 

21 

15 

15 

14405 


105 

33 

259 

87 

77 

63 

25 

4423 


91 

49 

505 

193 

125 

167 

127 

22916 


43 

29 

95 

67 

63 

57 

7 

77254 


49 

75 

319 

115 

103 

191 

109 

49143 

40 

27 

27 

159 

89 

87 

71 

45 

90282 


5 

5 

33 

37 

23 

25 

21 

8413 


55 

79 

705 

69 

65 

817 

361 

6355 


163 

167 

483 

277 

219 

191 

209 

39085 


5 

5 

7 

13 

11 

5 

5 

75880 


41 

41 

249 

165 

117 

35 

31 

38894 


3 

3 

89 

71 

51 

15 

7 

23092 


273 

187 

703 

209 

143 

337 

141 

38733 

42 

3 

3 

39 

83 

13 

15 

5 

155591 


31 

31 

87 

297 

69 

79 

85 

27015 


39 

35 

477 

67 

221 

251 

259 

17179 


85 

85 

477 

175 

97 

231 

181 

43292 


15 

33 

689 

259 

251 

123 

119 

44184 


11 

9 

273 

151 

79 

121 

27 

18039 


77 

75 

709 

271 

143 

299 

227 

31753 


27 

21 

211 

99 

147 

77 

35 

314516 

44 

7 

7 

671 

239 

161 

147 

33 

296542 


5 

5 

297 

73 

47 

5 

5 

88731 


59 

59 

273 

83 

87 

123 

73 

156877 


67 

41 

683 

79 

267 

255 

89 

161405 


43 

39 

395 

303 

173 

37 

53 

481254 


7 

13 

237 

235 

121 

105 

17 

572913 


11 

11 

125 

29 

149 

191 

159 

115119 


5 

5 

267 

105 

107 

125 

23 

126330 

46 

49 

29 

991 

447 

317 

159 

147 

140593 


25 

17 

167 

137 

105 

9 

17 

45621 


55 

51 

431 

137 

93 

127 

59 

467641 


51 

55 

731 

31 

13 

147 

89 

361410 


13 

9 

175 

115 

147 

15 

21 

848229 


33 

33 

123 

91 

97 

39 

17 

193586 


33 

29 

351 

187 

173 

45 

31 

102096 


5 

5 

273 

89 

33 

41 

25 

737685 

48 

23 

23 

67 

93 

63 

41 

11 

1151380 


109 

65 

181 

207 

181 

65 

39 

869929 


11 

11 

27 

47 

59 

55 

15 

1640100 


5 

5 

25 

129 

21 

29 

29 

256597 


47 

29 

389 

269 

117 

183 

31 

925482 


19 

19 

115 

67 

67 

11 

23 

264846 


31 

35 

689 

405 

227 

147 

99 

550320 


15 

15 

253 

65 

43 

91 

67 

769463 

mean 




110 

83 

71 

42 

92978 

wins 




0 

1 

5 

8 

0 
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Table 2. Time required to solve each instance 


size 

mostviol 

allviol 

allcst 

lp4 

lp8 

freq4 

freq8 

gurobi 

queries 

QAL 

38 

0.955 

0.939 

1.063 

1.96 

2.537 

1.627 

1.2 

0.849 

37 

-2 


1.881 

1.823 

1.721 

1.574 

0.799 

1.034 

1.46 

0.181 

81 

-8 


2.833 

2.987 

0.762 

1.919 

2.272 

2.432 

1.852 

0.362 

78 

-5 


0.56 

0.458 

0.728 

0.388 

0.456 

0.299 

0.458 

0.761 

40 

12 


2.182 

0.909 

3.078 

1.18 

1.451 

1.72 

1.321 

0.264 

186 

-3 


3.513 

2.529 

6.891 

2.934 

3.267 

3.415 

4.343 

1.68 

271 

-3 


2.262 

1.658 

0.898 

1.192 

1.728 

1.075 

0.342 

3.496 

23 

137 


3.108 

4.313 

4.055 

2.143 

2.704 

2.369 

2.896 

2.573 

319 

1 

40 

1.747 

1.72 

2.305 

1.581 

3.478 

2.205 

2.609 

6.381 

322 

15 


1.082 

1.099 

0.462 

0.786 

0.786 

0.475 

0.606 

0.53 

44 

2 


2.895 

3.108 

7.949 

3.032 

3.106 

11.263 

5.38 

0.465 

79 

-31 


9.201 

12.972 

7.575 

5.459 

5.87 

5.331 

6.912 

2.452 

370 

-8 


0.982 

1.088 

0.115 

0.12 

0.116 

0.166 

0.209 

3.991 

15 

258 


2.487 

2.534 

3.449 

2.332 

1.889 

0.821 

0.982 

2.413 

68 

23 


0.556 

0.57 

1.144 

1.714 

1.626 

0.584 

0.927 

1.444 

13 

68 


5.398 

3.613 

15.422 

4.691 

3.236 

7.188 

5.167 

3.055 

1121 

-0 

42 

0.716 

0.682 

0.39 

0.989 

0.399 

0.881 

0.323 

11.271 

17 

644 


4.602 

3.565 

3.322 

10.625 

2.315 

1.393 

2.127 

2.123 

104 

7 


5.464 

5.693 

13.078 

2.092 

8.149 

4.903 

4.339 

1.516 

147 

-4 


7.196 

7.144 

12.966 

4.138 

3.683 

3.428 

4.195 

3.102 

275 

-1 


2.78 

5.195 

10.686 

7.254 

9.431 

3.734 

4.188 

2.852 

57 

1 


4.248 

1.115 

2.981 

2.692 

2.19 

1.649 

1.68 

1.444 

36 

9 


4.948 

5.213 

9.818 

5.82 

5.182 

5.295 

4.964 

1.998 

231 

-13 


3.57 

2.296 

3.93 

2.276 

3.989 

2.701 

2.576 

19.461 

386 

45 

44 

2.63 

2.525 

9.813 

5.021 

3.768 

3.611 

1.99 

20.968 

79 

240 


1.423 

1.387 

4.59 

1.919 

1.607 

0.214 

0.275 

7.166 

18 

386 


3.246 

3.178 

6.486 

1.467 

2.295 

2.487 

1.273 

9.517 

117 

70 


3.756 

3.037 

22.551 

2.546 

5.31 

5.413 

5.486 

12.428 

223 

44 


5.699 

5.428 

5.237 

6.491 

5.841 

1.545 

3.787 

34.576 

149 

222 


1.89 

2.113 

3.566 

3.161 

2.728 

2.066 

1.654 

39.706 

56 

680 


1.906 

1.934 

4.101 

0.789 

5.77 

3.657 

3.338 

8.019 

57 

127 


1.591 

1.568 

4.824 

2.296 

3.413 

2.25 

2.047 

10.254 

21 

414 

46 

5.281 

3.116 

26.239 

12.163 

8.115 

3.593 

8.537 

13.92 

182 

59 


3.341 

2.543 

3.244 

3.006 

4.047 

0.404 

0.991 

4.293 

36 

108 


6.037 

5.757 

4.988 

3.342 

2.689 

2.102 

2.844 

33.412 

257 

122 


7.495 

7.725 

10.645 

0.497 

0.195 

5.247 

2.955 

30.772 

53 

577 


1.708 

1.855 

2.402 

3.173 

3.353 

0.614 

1.409 

70.359 

47 

1484 


3.653 

3.668 

1.944 

1.99 

2.373 

1.113 

1.143 

13.246 

65 

187 


5.767 

4.422 

5.828 

4.132 

4.435 

2.199 

2.359 

10.551 

324 

26 


1.741 

1.781 

3.203 

1.902 

0.783 

1.878 

1.653 

52.928 

74 

705 

48 

3.317 

3.286 

1.291 

3.149 

1.045 

2.698 

0.969 

81.257 

36 

2230 


14.462 

10.672 

9.82 

6.138 

9.82 

7.238 

5.852 

84.355 

293 

268 


2.464 

2.45 

0.551 

0.999 

1.638 

1.604 

1.019 

108.53 

79 

1367 


1.771 

1.944 

0.48 

5.574 

1.414 

1.52 

1.232 

19.495 

42 

453 


10.927 

4.585 

11.684 

10.371 

3.837 

14.762 

5.156 

96.16 

1187 

78 


4.197 

4.24 

1.202 

1.841 

2.871 

0.56 

1.337 

22.136 

21 

1027 


6.55 

4.952 

20.775 

18.857 

12.067 

8.16 

8.582 

46.879 

162 

259 


2.835 

2.758 

3.807 

1.065 

1.376 

2.738 

2.77 

48.871 

123 

389 

mean 

2.601 

2.617 

3.403 

2.437 

2.392 

1.968 

1.927 

6.605 

- 

- 

wins 

2 

4 

4 

6 

3 

11 

7 

11 

- 

- 
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